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I • In tr od uct ion 

The use of uniformly or non-unif ormly distributed 
pseudorandom numbers in systems simulation, statistical 
sampling experiments and analytical Monte Carlo work is by 
now well established. Numerous algorithms exist for 
producing such numbers from various distributions; for 
summaries of common techniques, see Knuth [5], Gaver and 
Thompson [2] or Ahrens and Dieter [1]. 

The user of pseudorandom numbers is usually not 
concerned with the details of the algorithm employed but 
rather with the results; a good algorithm, then, is one 
which is fast, uses minimum computer memory and produces 
numbers with satisfactory statistical properties. The 
search for statistically competent algorithms for 
pseudorandom numbers has resulted in the specif ication of 
many so-called "exact” generators, that is those whose 
deviation from the true distribution concerned is the result 
of computer rounding errors rather than any defect in the 
method itself. Such methods for nonuniform random numbers 
are often based on the assumption that "good" uniform 
numbers are available from an independent generator. 

Exact generators for nonuniform pseudorandom numbers 
are often quite complex and so assembly-level coding is 
often resorted to when implementing them in order to meet 
the computer time and memory constraints on a good 
algorithm. An example is the LLRANDOM package developed at 
the Naval Postgraduate School by G.P. Learraonth and P.A.W. 
Lewis and described in [7]; it produces pseudorandom numbers 
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from uniform, normal and exponential distributions. This 
report describes an extension to the LLRANDOM package for 
Cauchy and gamma distributed numbers. 



The C a u ch y distr ib ution has density function 
(1) f (x) = 1 1__, , - “ < x < ~ , 

TT T + X? 

and distribution function 



F(x) 



1 + 
2 



1 tan -1 x. 

TT 



While the shape of the Cauchy density resembles the normal 
density, the tails are much heavier; in fact, Cauchy random 
variables have no expectation and an infinite variance. The 
density has mode at zero and often in applications the 
variates are often shifted by a location parameter T or 
scaled by multiplying by a scale parameter S. Because of 
the heavy tails, Cauchy variates might find application ag a 
"pathological" case in a systems simulation study as well as 
in statistical sampling experiments for robust estimation 
techniques. See Chapter 16 of Johnson and Kotz [4] for 
further details on the Cauchy distribution. 



The gamma distribution with shape parameter A and scale 
parameter s has the density function 



A A- 1 -sx 

(2) f(x) = s x e , 

rOT — 



where r (A) is Euler's gamma function 



(3) 



r (A) 



A- 1 -x 
x e dx 



Note that r (n) = (n-1) ! when n is a non-negative integer. 

If the random variable X has density (2) then 



E[X] = A / s , 



2 



V[X] = A / s 2 



When A = 1, X has the exponential distribution while X, 
suitably scaled, has an asymptotically normal distribution 
as A -> oo . 

We note that if X has a r (A,1) distribution then X/s 
has a r <A,s) d istribution, so we may set s = 1 in (2) as 
far as the generating algorithm is concerned. The output 
from the generator may then be appropriately scaled. 

Gamma random variables are used in a wide variety of 
applications: for analytical modeling, in reliability 



theory and for 


statistical 


testing (the 


chi- 


squared random 


variable with 


n degrees 


of freedom 


has 


the T (n, 

2 




distribution) . 


See [ 6 ] or 


Chapter 17 


of 


[ 4 ] for 


more 



details. 
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II. U§e of the Subroutines 



This extension to LLRANDOM is composed of two 
independent IBM System/360 Assembler-coded subroutines: 

CAUCHY for Cauchy-distributed variates and GAMA for gamma 
variates. The name GAMA was chosen so. as not to conflict 
with the IBM mathematical library subprogram GAMMA which 
computes the gamma function (3) . 

The basic conventions for using GAMA and CAUCHY are the 
same as in the LLRANDOM package: the invoking statements 

CALL CAUCHY ( IX, X, N ) 
and CALL GAMA ( A, IX, X, N ) 

will result in a vector X(1), ... , X(N) of Cauchy or 

r (A, 1.0) pseudorandom variates, respectively. The arg- 
ument IX is, in both cases, an integer seed to be used in 
the multiplicative congruential uniform generator employed 
by LLRANDOM. IX should be initialized just once in the 
calling program to some positive integer value and should 
not be altered thereafter. 

The subroutine GAMA requires a source for normal and 
exponential deviates; these are obtained directly from the 
LLRANDOM package and so the statement "CALL OVFLOW" must 

appear once in the calling program to initialize LLRANDOM. 

As mentioned previously, the output from GAMA must be scaled 
if the scale parameter is other than one; the following set 
of statements will thus be required to generate a vector of 
100 chi-squared variates with seven degrees of freedom: 

DIMENSION X (100) 

CALL OVFLOW 

IX = 13726 

• • • 

CALL GAMA ( 3.5, IX, X, 100 ) 
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DC 50 I = 1,100 
X(I) = 2.0 * X(I) 
50 CONTINUE 

• • • 

END 



Cauchy variates are also often modified by location and 
scale parameters; since no expectations exist, however, we 
cannot refer to these parameters in terms of mean or 
variance. Subroutine CAUCHY is completely independent of 
LLRANDOM or any other subroutines so that the "CALL OVFLOW" 
statement is not necessary in this case. To use CAUCHY to 
produce a single variate C with location parameter T and 
scale parameter S we may use the statements 

• • • 

IX = 217663541 

CALL CAUCHY ( IX, C, 1 ) 

C = S * C + T 

END 



Just as in LLRANDOM, linkage overhead between the 
calling program and GAMA or CAUCHY will be minimized if a 
vector of several variates is obtained at the same time 
instead of just a single one. The gain in this case can be 
as much as 50 microseconds per variate in average generation 
time, an improvement of up to 50%. In GAMA, several 
constants must be calculated for each different value of the 
shape parameter A; these constants are saved between calls 
so that they need not be recomputed. It will thus be more 
efficient to get several gamma variates with the same shape 
parameter before changing the A value, especially when A > 
3.0 when the setup computations are extensive (see lines 
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174-246 of the program listing) . 



Note that the techniques used in GAMA and CAUCHY make 
use of so-called rejection methods so that the number of 
uniform (or exponential or normal) deviates needed to 
generate a single output deviate is random. When normal or 
exponential deviates are required by GAMA from LL RANDOM a 
vector of 10 deviates is called for; since not all of these 
may be used at the time they are generated, the balance are 
saved for the next call to GAMA. Thus, reinitializing the 
seed IX to its original value will not in general result in 
an exact repetition of the generated gamma sequence since 
the first few deviates will use the old normal or 
exponential deviates from the previous sequence. To achieve 
an exact repetition, the generator must be forced to repeat 
the initialization computations for the desired A value; at 
this time any remaining variates from LLRANDOM are 
discarded. An example of this might be 

DIMENSION G (100) 

CALL OVFLOW 
IX = 12345 



CALL GAMA ( A, IX, G, 100 ) 

• • • 

C REINITIALIZE GAMMA SEQUENCE 

CALL GAMA ( 1.0, IX, G, 1 ) 
IX = 12345 



CALL GAMA ( A, IX, G, 100 ) 

• • • 

END 



CAUCHY requires 552 bytes and, as mentioned previously, 
is completely independent of any other subprograms. CAUCHY 
uses the LLRANDOM multiplicative congruential uniform 
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generator but. this is coded in line when needed so as to 
preserve CAUCHY’s independence. The average generation time 
per variate for subroutine CAUCHY on a System/360 Model 67 
under OS/MVT was 67.5 microseconds when variates were 
generated in vectors of 100. The generation of variates one 
at a time increased the average time to 119.3 microseconds 
per variate. 

Subroutine GAMA itself uses only 1988 bytes of memory 
but since it calls on LLRANDOM the total core requirement is 
9342 bytes: 



GAMA 1988 bytes 
LLRANDOM 6189 bytes 
Required IBM Functions 2265 bytes 
Total 9342 bytes 

Timing the gamma generator on a System/360 Model 67 was 
carried out using the TIME macro; Table 1 summarizes the 
observed times as a function of the shape parameter, A. 
Note that since special methods are employed when A is 0.5, 
1.0, 1.5, 2.0 or 3.0, the times in these cases are 
considerably shorter than times for nearby values of A. 
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Shape Parameter Algorithm Vector of 100 Single 



A 




Variates 


Variate 


0 . 1 


GS 


324.0 


364.0 


0.3 


GS 


367.0 


402.5 


0.5 


GA 


70.4 


207.7 


00 

• 

o 


GS 


439.8 


551.2 


0.9 


GS 


459.0 


611.0 


1.0 


GA 


68.7 


158.9 


1 . 2 


GF 


300.1 


385.0 


1 . 4 


GF 


306. 1 


441 .0 


1.5 


GA 


141.7 


215.8 


1.8 


GF 


343.6 


390.8 


2.0 


GA 


142.5 


203.6 


2. 1 


GF 


396.1 


450.8 


2.5 


GF 


434.7 


468.5 


2.9 


GF 


444.5 


496.6 


3.0 


GA 


206.7 


237.1 


3. 1 


GO 


341.5 


435.8 


3.5 


GO 


336.2 


373.4 


4.0 


GO 


332.4 


420.7 


5.0 


GO 


307.7 


363.2 


8.0 


GO 


293. 1 


371.3 


10.0 


GO 


289.4 


312.5 


20.0 


GO 


238.2 


321.6 


50.0 


GO 


197.7 


284.2 


100. 0 


GO 


178.4 


220.0 


1000.0 


GO 


166.7 


177.0 


10000.0 


GO 


136.4 


169.8 


100000.0 


GO 


152.5 


235.8 



Table 1. Average generation times (microseconds) for gamma 
variates using subroutine GAMA. 
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III. Desc ri pt ion of the Algo ri thms 



This section describes the actual algorithms used in 
CAUCHY and GAMA. An understanding of the algorithms is not 
necessary for use of the package but they are set forth here 
both in the interest of completeness and in an effort to 
document the programs more fully. A single algorithm 
suffices for the Cauchy generator while GAMA uses one of 
four algorithms, depending on the value of A. 

In the descriptions which follow, the letters U, N and 
E (with or without affixes) represent uniform, standard 
normal and unit exponential pseudorandom deviates, 
respectively. The phrase "Generate U" implies that U is the 
next sequential uniform variate in the linear congruential 
sequence; these variates are generated as needed by using 
the same multiplicative congruential scheme as used in 
LLRANDOM. The phrases "Generate N" or "Generate E" imply 
that normal or exponential variates are to be obtained by 
linking directly to LLRANDOM. 

A . Cauchy Genera tor 



The Cauchy generator is a combination decora position- 
rejection method (see Knuth [5]). The Cauchy density is 
decomposed, as in Figure 1, into three subdensities: a 

uniform density between 0 and 1 (f^), a wedge-shaped density 

(f^) and a long tailed density (f^) . 



The 
1/ it ; in 
density 



uniform density f^ is sampled with probability 
this case a uniform (0, 1) variate is returned. The 
f is dealt with by using Marsaglia's almost-linear 
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density algorithm, just as in Knuth's Algorithm L [5]. The 
density f is sampled with probability 1/2 - 1 /tt . The tail 
density f is sampled by a rejection method with probability 
1/2. The majorizing density for f^ is g(x) = 1 / x 2 , which 
is the density of the reciprocal of a uniform (0,1) variate. 

Algorithm C below uses the fact that in the prime 

modulus congruential random number generator used in 

LLRANDOH the low order bits are uniformly distributed so 

that b and b select the proper sub-distribution in Step 1. 
12 

This will not in general be the case for other congruential 
pseudo-random number generators. 




Figure 1. Decomposition of the Cauchy Density Function. 
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Algorith m C. Cauchy variates. 

1. (Select subdensity) Generate U, setting aside the two 

low order bits b and b . If b = 1, go to Step 6. 

12 1 

2. (Sample box) If II < 0.6366197724 = 2/ v , generate a 

* * 

new variate U , set x = U and go to Step 8. 

3. (Sample wedge) Generate new variates U and U . If 0 

1 2 1 

> U , exchange U and U . Set x = U . 

2 12 1 

4. (Easy rejection) If < 0.8284271247 =2/2-2, go 
to Step 8. 

5. (Hard rejection) If (J - U < 1 - x 2 ( 2/2 - 2) , go to 

2 1 T““x 7 

Step 8, otherwise go back to Step 3. 

6. (Sample tail) Set x = 1 / U. 

if if 

7. (Tail rejection) Generate a new variate U . If U < 

_x£_ go to Step 8, otherwise generate a new U and go 
T + x 7 

back to Step 6. 

8. (Random sign) If b = 1 set x = - x. Deliver x as the 
generated deviate. 



It should be noted that there are several other methods 
for generating Cauchy variates: the ratio of independent 

standard normal deviates has the Cauchy distribution, as 
does the quantity 

X = tan [ tt (U - 1) ], 

2 

where U is unif orm (0 , 1 ) . These methods are both substant- 
ially slower than algorithm C, but another new method has an 



1 1 



average time comparable to Algorithm C and is much easier to 
program. This second method requires an average of 2.55 
uniform random variates per Cauchy variate (as compared with 
2.47 for algorithm C) and it needs about 69 microseconds per 
variate on the System/360 Model 67. It is possible, 
however, that Algorithm CR will be better than algorithm C 
in some other implementation. 

The method is essentially the technique devised by von 
Neumann to generate a random variate sin (J, where U is 
uniform between 0 and 2n . Such variates are used in the 
polar method for generating normal random variables [8]. It 
does not seem to have been recognized that the method also 
generates tan U, which is the required Cauchy variate. 



Alg orit hm CR. Cauchy variates, ratio method. 

1. (Get uniforms) Generate U and U . Set Y = 2 U 

12 11 



2 . 



3. 



and Y = 2 U - 1. 



2 2 

(Rejection test) If Y* + y^ > 1 go back to Step 

(Take ratio) Deliver x = Y / Y . 

12 



1 . 



1 



B. Gamma Ge ner at or GS : A < .K 0 



This method is due to Ahrens and is set forth in [ 1 ]. 
It is applicable only to values of A less than one and is 
markedly superior in execution time to the method of Johnk 
[3], which is the usual technique for generating variates of 
this type. 

The method is a rejection method employing two 
different tests, one of which is chosen at random for any 
given variate: the power transform of a uniform (0,1) 
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variate, U , is tested in the region 0 < x < 1, while a 

suitable exponential, E, is tested when x > 1. The 
advantage of this method lies in the limited use of the 
library subprograms for the exponential and logarithm; 
average times range from 300 to 400 microseconds as compared 
with 600 to 800 for Johnk's method. Further discussion and 
proofs may be found in [1]. 



Algorithm GS. Gamma variates, A < 1.0. 

1. (Select rejection test) Generate U and generate E and 

set P = e + A U. (Note that "e" is the base of the 
e 

natural logarithms.) If P < 1 go to Step 2, otherwise 
go to Step 3. 

1/A 

2. (Small x test) Set x = P . If x < E, deliver x, 

otherwise go back to Step 1. 

3. (Large x test) Set x = - In [ 1 { e_+_A - P } ]. If 

~K e 

( 1 - A) In x < E, deliver x, otherwise go back to Step 

1 . 

C . G^mma Genera tor GF: 1 . 0 < A < 3 .0 

A thus-far unpublished method devised by Professor G.S. 
Fishman of North Carolina University was communicated to the 
authors in private correspondence. It is valid for any A > 
1.0 but its efficiency in terms of average time goes down as 

v/S so it is applied in GAMA only in the range where it is 

superior to the Dieter- Ahrens method GO described below. 
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The method is a rejection method based on the following 
theorem. 

Theor em Let U be a uniform (0,1) random variable and let E 
be an exponential random variable with mean A. Let 

A-1 -X (1—1/A) - ( A — 1 ) 

= [ i 3 e 

If g (E) > U, then E has conditionally the gamma distribution 
with shape parameter A, i.e. 

A*" 1 — x 

f ( x | U < g (E) ) = x e 

e *TW 



Proof : 



x/A 

Unconditionally, E has density h (x) = 1 e 

Therefore, 

(4) f E (x|U<g(E) ) = hjxl_Pr jU<g.jE). p= xl . 

Now since 0 is uniformly distributed, 

Pr {U<g (E) | E=x} = g (x) 

as long as 0 < g (x) <1; that this is true for every x > 0 
may be readily verified by elementary calculus. Therefore, 

(5) Pr {U<g (E) } = E[ Pr{U<g(E)|E} ] 

= j g(x) h(x) dx 
A-1 -A 

= r (A) e A 
= C (A) 

Thus, in view of (4) , 
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f (x | U <g (E) 
E 



A- 1 -x 
x e 

~tw 



The efficiency of the generator is governed by the 
probability that a given variate will pass the rejection 
test, U < g (E) ; from (5) it will be seen that this 
probability is just C(A) . When A is large we have from 

Stirling's approximation that C(A) = /~2jT , so that the 

y I e* 

method oecomes more inefficient with increasing A, as noted 
above. 

A slight modification to the method suggested by the 
theorem improves the efficiency slightly and we obtain 



Ai 3 °Eithm GF. Gamma variates, 1.0 < A < 3.0. 

1. (Generate exponentials) Generate two 

exponential variates, E and E . 

12 

2. (Rejection test) If E < (A— 1 ) (E - In E 

2 1 1 

go back to Step 1. 

3. (Acceptance) Deliver x = A E^. 



independent 



1) then 



D. Gam ma Generato r GO; A > 3^0 

This method was originally developed by Dieter and 
Ahrens and is fully described in [ 1 ] together with several 
other gamma generation techniques. Algorithm GO does not 
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suffer the usual drawback of growing less efficient in 
generation time with increasing A; in fact, the method is 
more efficient for larger A values. 

The basic idea here is to take advantage of the 
asymptotic normality of the gamma distribution by doing most 
of the sampling from a normal distribution; the right hand 
tail is sampled, when necessary, using a rejection method 
with the exponential distribution. The method can be 
applied to values of A greater than 2.533, but it is not as 
efficient as Fishman’s technique for A < 3.0. 

As mentioned previously, this algorithm requires the 
computation of several constants which depend only on A and 
which may be saved between calls; these calculations are 
described in step 0 of the specification below. Further 
discussion, illustrations and proofs are given in [ 1 ]; the 
version of GO here differs in a few minor details from the 
original Dieter and Ahrens technique. 



Algorithm GO. Gamma variates, a > 3.0. 

0. (Calculate constants) Compute: 

m = A - 1 ; 



s 2 = 



/I s * 



A; 



s = /s J ; 



d = /Fs 7 ; b = d + m ; 

w = s 2 / m - 1; v = 2s 2 / (m /I ) ; 

c=b+lnsd- 2m- 3.7203285. 

E 



1 . (Select normal/exponential) Generate U. If U < 

0.0095722652 go to Step 7. 

2. (Normal sampling) Generate N and set x = sN + m. 

3. (Check trial value) If x < 0 or x > b go back to Step 2, 
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N 2 / 2. 



otherwise generate a new variate U and set S = 

If N > 0 go to Step 5. 

4. (Left-hand rejection) If 0 < 1 + S (vll - w) go to Step 
9, otherwise go to Step 6. 

5. (Right-hand rejection) If U < 1 - wS go to Step 9. 

6. (Final normal rejection) If In U < m In x + m - x + S 

m 

go to Step 9; otherwise go back to step 1. 

7. (Exponential) Generate E and E and set x = b(1 + E /d) . 

1 2 1 

8. (Exponential rejection) If m ( x - In x ) + c > E go 

15 m 2 

back to Step 1 . 

9. (End) Deliver x as the gamma variate. 

E. Ad Hoc Gamma Generators 

This set of algorithms is based on the well-known fact 
that the sura of independent gamma variates with shape 

parameters A and A and equal scale parameters has the 

1 2 

gamma distribution with shape parameter A + A and scale 

1 2 

parameter equal to that of the summands. We may thus 

generate a gamma variate with integer shape parameter K by 
taking the sum of K independent exponentials. This will be 
more efficient than the previously discussed methods 
(Algorithms GF and GO) for moderate values of K; for the 
System/360 we take K < 3 to apply this ad hoc technique. 

An obvious extension to this method is to allow for 
half-integral values of A by making use of the fact that the 
square of a standard normal random variable has the 
chi-squared distribution with one degree of freedom, i.e. 
N 2 /2 has the gamma distribution with unit scale parameter 
and A = 0.5. We use this extension for A = 0.5 or 1.5. 
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The resulting algorithm is then 



Algorithm GA. Gamma variates, integral or half-integral 

shape parameter A. 

1. (Find K) Set K = [A], where [A] denotes the integral 

part of A. Set X = 0 . If A - K = 0.5 set L = 1; if A - 

K = 0.0 set L = 0; otherwise Stop. (If the algorithm 

stops, an incorrect A value has been used.) 

2. (Generate exponentials) If K = 0 go to Step 3, 

otherwise generate K exponentials E , ... ,E and set 

1 K 

X = E + . . . ♦ E • 

1 K 

3. (Generate normal) If L = 0 go to Step 4 otherwise 

generate N and set X = X + N 2 /2. 

4. (Deliver X) X is the desired variate. 
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IV. Sum n ar y and Com ment s 



This work provides a convenient and useful extension to 
the LLRANDOM package, especially for users interested in 
statistical and reliability theory applications of digital 
simulation. The combination of the most efficient known 
gamma generation techniques with the new Cauchy method gives 
exceptionally good time characteristics at some cost in 
computer memory utilization. 

The work may be extended at once to the generation of 
several other types of random variables. For example, the 
beta distribution with parameters A and B may be sampled by 
taking gamma variates X^ and X^ with respective shape 

parameters A and B and delivering 

Z = X / ( X + X ) 

1 1 2 

as a beta variate. In this case considerable overhead in 

GAMA can result from shifting the shape parameter back and 

forth between A and B; for this reason obtaining vectors of 

gamma variates X and X is recommended, as in the following 

12 

example : 

DIMENSION XI (50), X2 (50) , Z (50) 

CALL GAMA ( A, IX, XI, 50 ) 

CALL GAMA ( B, IX, X2, 50 ) 

DO 405 I = 1,50 

Z(I) = XI (I) / ( XI (I) + X 2 (I) ) 

405 CONTINUE 

END 
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The t- Pis tribat ion may be sampled as the ratio of a 
standard normal and an independent chi-sguared random 
variate, while the F -Di stribu tio n may be obtained by taking 
the ratio of two independent chi-sguared variates divided by 
their respective degrees of freedom. (See pages 4 and 5 for 
an example of the generation of chi-sguared variates.) 
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**** GAMMA deviate generator **** 
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** GAMMA DEVIATE GENERATOR **** 
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*❖*4= gamma deviate generator 
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**** GAMMA DEVIATE GENERATOR **** 
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**** GAMMA DEVIATE GENERATOR 
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